May, 2018
if(!require(OpenStreetMap)){install.packages("OpenStreetMap");
library(OpenStreetMap)}
## Loading required package: OpenStreetMap
if(!require(ggplot2)){install.packages("ggplot2"); library(ggplot2)}
## Loading required package: ggplot2
* 웹에서 제공하는 지도정보를 사용할 수 있는 패키지 * `OpenStreetMap`을 이용하기 위한 패키지 호출
map = OpenStreetMap::openmap(upperLeft = c(43, 119),
lowerRight = c(33, 134), type = 'bing')
plot(map)
* upperLeft에는 지도 왼쪽 위 모서리의 위경도 좌표를 입력
* lowerRight에는 지도 오른쪽 아래 모서리의 위경도 좌표를 입력
* ggplot2 패키지의 autoplot 함수로 openmap 을 그릴 수 있음
nm = c("osm", "mapbox", "stamen-toner",
"stamen-watercolor", "esri", "esri-topo",
"nps", "apple-iphoto", "osm-public-transport")
par(mfrow=c(3,3), mar=c(0, 0, 0, 0), oma=c(0, 0, 0, 0))
for(i in 1:length(nm)){
map <- openmap(c(43,119),
c(33,134),
minNumTiles = 3,
type = nm[i])
plot(map, xlab = paste(nm[i]))
}
par(mfrow = c(1, 1))
* openmap 함수의 type의 종류는 bing 외에도 많이 있다.
+ nm 벡터에 웹지도의 소스에 대한 이름이 들어 있다.
map1 <- openmap(c(43.46,119.94),
c(33.22,133.98))
plot(map1)
abline(h = 38, col = 'blue')
str(map1)
## List of 2 ## $ tiles:List of 1 ## ..$ :List of 5 ## .. ..$ colorData : chr [1:378378] "#FAF9F6" "#585857" "#FAF9F6" "#FAF9F6" ... ## .. ..$ bbox :List of 2 ## .. .. ..$ p1: num [1:2] 13351660 5382253 ## .. .. ..$ p2: num [1:2] 14914585 3924542 ## .. ..$ projection:Formal class 'CRS' [package "sp"] with 1 slot ## .. .. .. ..@ projargs: chr "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs" ## .. ..$ xres : int 594 ## .. ..$ yres : int 637 ## .. ..- attr(*, "class")= chr "osmtile" ## $ bbox :List of 2 ## ..$ p1: num [1:2] 13351660 5382253 ## ..$ p2: num [1:2] 14914585 3924542 ## - attr(*, "zoom")= int 6 ## - attr(*, "class")= chr "OpenStreetMap"
map1$tiles[[1]]$projection
## CRS arguments: ## +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 ## +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs
if(!require(sp)){install.packages("sp"); library(sp)}
map_p <- openproj(map1, projection = CRS("+proj=longlat"))
str(map_p)
## List of 2 ## $ tiles:List of 1 ## ..$ :List of 5 ## .. ..$ colorData : chr [1:393336] NA NA NA NA ... ## .. ..$ bbox :List of 2 ## .. .. ..$ p1: num [1:2] 119.8 43.6 ## .. .. ..$ p2: num [1:2] 134.1 33.1 ## .. ..$ projection:Formal class 'CRS' [package "sp"] with 1 slot ## .. .. .. ..@ projargs: chr "+proj=longlat +ellps=WGS84" ## .. ..$ xres : int 607 ## .. ..$ yres : int 648 ## .. ..- attr(*, "class")= chr "osmtile" ## $ bbox :List of 2 ## ..$ p1: num [1:2] 119.8 43.6 ## ..$ p2: num [1:2] 134.1 33.1 ## - attr(*, "zoom")= int 6 ## - attr(*, "class")= chr "OpenStreetMap"
plot(map_p) abline(h = 38, col = 'blue')
map_p <- openproj(map1, projection =
CRS("+proj=utm +zone=52N + datum=WGS84"))
plot(map_p)
abline(h = 38, col = 'blue')
if(!require(sp)){install.packages("sp"); library(sp)}
a <-data.frame(lon = seq(100,140,by = 0.1),
lat = 38)
sp::coordinates(a) = ~ lon + lat
str(a)
## Formal class 'SpatialPoints' [package "sp"] with 3 slots ## ..@ coords : num [1:401, 1:2] 100 100 100 100 100 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:401] "1" "2" "3" "4" ... ## .. .. ..$ : chr [1:2] "lon" "lat" ## ..@ bbox : num [1:2, 1:2] 100 38 140 38 ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:2] "lon" "lat" ## .. .. ..$ : chr [1:2] "min" "max" ## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot ## .. .. ..@ projargs: chr NA
sp::proj4string(a) = "+proj=longlat"
#a@proj4string = CRS("+proj=longlat")
str(a)
## Formal class 'SpatialPoints' [package "sp"] with 3 slots ## ..@ coords : num [1:401, 1:2] 100 100 100 100 100 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:401] "1" "2" "3" "4" ... ## .. .. ..$ : chr [1:2] "lon" "lat" ## ..@ bbox : num [1:2, 1:2] 100 38 140 38 ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:2] "lon" "lat" ## .. .. ..$ : chr [1:2] "min" "max" ## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot ## .. .. ..@ projargs: chr "+proj=longlat +ellps=WGS84"
a_tf = spTransform(a, CRS("+proj=utm +zone=52N + datum=WGS84"))
str(a_tf)
## Formal class 'SpatialPoints' [package "sp"] with 3 slots ## ..@ coords : num [1:401, 1:2] -2069280 -2060288 -2051296 -2042306 -2033316 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:401] "1" "2" "3" "4" ... ## .. .. ..$ : chr [1:2] "lon" "lat" ## ..@ bbox : num [1:2, 1:2] -2069280 4205815 1467197 4626502 ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:2] "lon" "lat" ## .. .. ..$ : chr [1:2] "min" "max" ## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot ## .. .. ..@ projargs: chr "+proj=utm +zone=52N +ellps=WGS84"
plot(map_p) points(a_tf@coords[,1], a_tf@coords[,2], type = 'l', col = 'blue')
* `openmap`으로 그린 지도위에 새로운 시각정보를 그리기 위해서
`openmap`이 생성한 좌표계를 알아야 한다.
+ 위경도
+ 구면좌표의 지도위로의 정사형 방법
* 새로운 시각정보를 표시하는 위치는 이미 그려진 지도의 좌표계를 따른다.
* 좌표계를 변환하는 라이브러리인 `sp`라이브러리를 사용하자.
if(!require(sp)){install.packages("sp"); library(sp)}
if(!require(mapplots)){install.packages("mapplots"); library(mapplots)}
map = openmap(upperLeft = c(43, 119),lowerRight = c(33, 134))
seoul_loc = geocode('Seoul')
coordinates(seoul_loc) = ~lon + lat
proj4string(seoul_loc) = "+proj=longlat +datum=WGS84"
seoul_loc_Tf = spTransform(seoul_loc,
CRS(as.character(map$tiles[[1]]$projection)))
plot(map)
add.pie(z=1:2,labels = c('a','b'),
x = seoul_loc_Tf@coords[1],
y = seoul_loc_Tf@coords[2], radius = 100000)
* `geocode` 함수를 이용하여 seoul의 위경도 좌표 호출
* `coordinates(seoul_loc) = ~lon + lat`를 이용하여 `sp`라이브러리에서 처리가능한
데이터 형으로 변환
* `sp` 패키지의 `spTransform()` 을 이용, openmap의 좌표계로 위경도 좌표를 변환
+ `map$tiles[[1]]$projection`에 이미 읽어들인 openStreetmap의 좌표정보가 있음.
+ `CRS` 함수를 이용하여 `CRSobj` 데이터 형으로 변환
if(!require(ggmap)){install.packages("ggmap"); library(ggmap)}
* ggmap을 이용하기 위한 패키지 호출 * 휴스턴의 범죄 데이터의 시각화 작업
data(crime) head(crime, 2)
## time date hour premise offense beat block ## 82729 2010-01-01 15:00:00 1/1/2010 0 18A murder 15E30 9600-9699 ## 82730 2010-01-01 15:00:00 1/1/2010 0 13R robbery 13D10 4700-4799 ## street type suffix number month day location ## 82729 marlive ln - 1 january friday apartment parking lot ## 82730 telephone rd - 1 january friday road / street / sidewalk ## address lon lat ## 82729 9650 marlive ln -95.43739 29.67790 ## 82730 4750 telephone rd -95.29888 29.69171
* crime 데이터에는 범죄 발생 시기, 위치, 종류 등이 저장
violent_crimes = subset(crime,
offense != "auto theft" &
offense != "theft" &
offense != "burglary")
violent_crimes$offense <- factor(violent_crimes$offense,
levels = c("robbery", "aggravated assault", "rape", "murder"))
violent_crimes = subset(violent_crimes,
-95.39681 <= lon & lon <= -95.34188 &
29.73631 <= lat & lat <= 29.78400)
* 범죄 데이터를 필요한 것만 선택해서 네가지의 범죄 종류를 factor 형태로 저장
+ `levels`:범죄의 중대성 순서대로 level을 정했음
* Houston의 downtown에서 발생한 범죄에 한정
HoustonMap = qmap("houston", zoom = 14,
color = "bw", legend = "topleft")
HoustonMap + geom_point(aes(x = lon, y = lat,
colour = offense, size = offense),
data = violent_crimes)
* `ggmap` 패키지의 `qmap` 함수를 이용하여 Houston의 지도 저장
+ `zoom`: 확대정도
+ `color`: 지도의 종류
* `geom_point`를 이용해서 저장된 지도 위에 `violent_crimes`를 점으로 시각화
+ `colour = offense`: 색깔을 이용해서 범죄의 종류를 구분
+ `colour = offense`: 범죄의 종류를 점의 크기에 적용
중대한 범죄의 레벨값을 높게 정했음
* Houston 지도 위의 범죄 데이터를 점으로 표기 <br>
* geom_density2d 함수를 이용하여 범죄 밀도 시각화
+ `size`: 선의 굵기
+ `bins`: 축방향의 bin 개수
HoustonMap +
geom_point(aes(x = lon, y = lat,
colour = offense, size = offense),
data = violent_crimes) +
geom_density2d(aes(x = lon, y = lat), size = 0.2 , bins = 4,
data = violent_crimes)
## Crimes map
* `stat_density2d` 을 이용한 density plot의 색깔 칠하기
load('airport.Rdata')
head(airport_krjp)
## city country iata lat lon Freq ## 1 Tokyo Japan NRT 35.76472 140.3864 100 ## 2 Matsumoto Japan MMJ 36.16676 137.9227 7 ## 3 Ibaraki Japan IBR 36.18108 140.4154 4 ## 4 Iwojima Japan IWO 24.78400 141.3227 NA ## 5 Nanki-shirahama Japan SHM 33.66222 135.3644 2 ## 6 Obihiro Japan OBO 42.73333 143.2172 6
* airport.Rdata는 `airport_krjp`, `link_krjp` 데이터를 포함 * airport_krjp는 공항의 데이터
head(link_krjp)
## iata lat lon group ## 1 CJJ 36.71660 127.4991 1 ## 2 CJU 33.51131 126.4930 2 ## 3 CJU 33.51131 126.4930 3 ## 4 CJU 33.51131 126.4930 4 ## 5 FUK 33.58594 130.4507 5 ## 6 GMP 37.55831 126.7906 6
* link_krjp는 비행의 출발지와 도착지의 데이터 * link_krjp 데이터에서 하나의 비행경로는 하나의 group으로 표기
map = ggmap(get_googlemap(center = c(lon=134, lat=36),
zoom = 5, maptype='roadmap', color='bw'))
map + geom_line(data=link_krjp,aes(x=lon,y=lat,group=group),
col='grey10',alpha=0.05) +
geom_point(data=airport_krjp[complete.cases(airport_krjp),],
aes(x=lon,y=lat, size=Freq), colour='black',alpha=0.5) +
scale_size(range=c(0,15))
* google map에서 지도를 가져와서 map 변수에 저장
+ `google map`은 기본적으로 위경도 좌표계를 사용하므로 좌표변환이 필요없음
* 지도에 `geom_line` 함수를 이용하여 비행 경로를 시각화
* `geom_point` 함수를 이용하여 점의 크기에 따라 공항 이용 빈도를 시각화
* 지도에 몇 개의 점에서 관측치를 표시
* 관측되지 않은 지점의 예측값을 표시하고 싶은 경우
+ 공간적 내삽 (spatial interpolation)
* Kriging 을 이용한 서울시 미세먼지 농도의 시각화 작업
if (!require(sp)) {install.packages('sp'); library(sp)}
if (!require(gstat)) {install.packages('gstat'); library(gstat)}
if (!require(automap)) {install.packages('automap'); library(automap)}
if (!require(rgdal)) {install.packages('rgdal'); library(rgdal)}
if (!require(e1071)) {install.packages('e1071'); library(e1071)}
if (!require(dplyr)) {install.packages('dplyr'); library(dplyr)}
if (!require(lattice)) {install.packages('lattice'); library(lattice)}
if (!require(viridis)) {install.packages('viridis'); library(viridis)}
데이터 정제 및 시각화를 위한 library 호출
seoul032823 <- read.csv ("seoul032823.csv")
head(seoul032823)
## X.1 X ID time PM10 LAT LON day ## 1 671 671 111121 2012032823 127 37.56464 126.9760 20120328 ## 2 2879 2879 111123 2012032823 121 37.57203 127.0050 20120328 ## 3 5087 5087 111131 2012032823 127 37.54031 127.0051 20120328 ## 4 7295 7295 111141 2012032823 143 37.54464 127.0957 20120328 ## 5 9503 9503 111142 2012032823 140 37.54311 127.0411 20120328 ## 6 11711 11711 111151 2012032823 128 37.58495 127.0943 20120328
* 서울시에서 제공하는 미세먼지 데이터 읽기
skorea <- raster::getData(name ="GADM", country= "KOR", level=2)
# skorea <- readRDS("KOR_adm2.rds")
head(skorea,2)
## class : SpatialPolygonsDataFrame ## features : 2 ## extent : 128.9957, 129.0844, 35.12158, 35.28471 (xmin, xmax, ymin, ymax) ## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 ## variables : 15 ## names : OBJECTID, ID_0, ISO, NAME_0, ID_1, NAME_1, ID_2, NAME_2, HASC_2, CCN_2, CCA_2, TYPE_2, ENGTYPE_2, NL_NAME_2, VARNAME_2 ## min values : 1, 213, KOR, South Korea, 1, Busan, 1, Buk, , NA, , Gu, District, 부산진구| 釜山鎭區, ## max values : 2, 213, KOR, South Korea, 1, Busan, 2, Busanjin, , NA, , Gu, District, 북구| 北區,
class(skorea)
## [1] "SpatialPolygonsDataFrame" ## attr(,"package") ## [1] "sp"
head(skorea@polygons[[1]]@Polygons[[1]]@coords, 3)
## [,1] [,2] ## [1,] 128.9957 35.17683 ## [2,] 128.9983 35.19250 ## [3,] 129.0003 35.20673
* GADM 데이터 구조 파악
+데이터의 구조가 'sp' 패키지에서 제공하는 S4-class 구조임을 확인
+head 함수로 확인할 수 없었던 좌표계를 @ 을 이용하여 확인
if (!require(broom)) {install.packages('broom'); library(broom)}
## Loading required package: broom
skorea <- broom::tidy(skorea)
## Regions defined for each Polygons
class(skorea)
## [1] "data.frame"
head(skorea,3)
## long lat order hole piece group id ## 1 128.9957 35.17683 1 FALSE 1 1.1 1 ## 2 128.9983 35.19250 2 FALSE 1 1.1 1 ## 3 129.0003 35.20673 3 FALSE 1 1.1 1
* 데이터의 구조를 data.frame 형태로 변경
ggplot() + geom_map(data= skorea, map= skorea,
aes(map_id=id,group=group),fill=NA, colour="black") +
geom_point(data=seoul032823, aes(LON, LAT, col = PM10),alpha=0.7) +
labs(title= "PM10 Concentration in Seoul Area at South Korea",
x="Longitude", y= "Latitude", size="PM10(microgm/m3)")
* geom_map 함수를 이용하여, GADM 데이터를 지도위에 plotting<br> * geom_point 함수를 이용하여, 미세먼지 데이터를 점으로 표기<br> * 보다 나은 시각화를 위해 점과 점 사이에 미세먼지 값의 내삽
class(seoul032823)
## [1] "data.frame"
coordinates(seoul032823) <- ~LON+LAT class(seoul032823)
## [1] "SpatialPointsDataFrame" ## attr(,"package") ## [1] "sp"
* 데이터 구조 변경
+Kriging 방법 적요을 위해 데이터를 'sp dataframe' 형태로 변경
+ 위경도 정보는 ``seoul032823@coords``에 저장되어 있음
LON.range <- c(126.5, 127.5) LAT.range <- c(37, 38) seoul032823.grid <- expand.grid( LON = seq(from = LON.range[1], to = LON.range[2], by = 0.01), LAT = seq(from = LAT.range[1], to = LAT.range[2], by = 0.01))
* Kriging을 위한 grid 생성
+ 새로운 좌표계의 범위를 계산
+ 좌표계의 범위를 이용하여 grid 생성
plot(seoul032823.grid) points(seoul032823, pch= 16,col="red")
* 참고: grid 위에서 데이터 시각화
coordinates(seoul032823.grid)<- ~LON+LAT ## sp class gridded(seoul032823.grid)<- T plot(seoul032823.grid) points(seoul032823, pch= 16,col="red")
* grid 데이터의 구조 변경 및 격자로 변경후 시각화
seoul032823_OK <- autoKrige(formula = PM10~1,
input_data = seoul032823,
new_data = seoul032823.grid )
## [using ordinary kriging]
* `autoKrige` 함수를 통해 만들어직 격자데이터 `seoul032823.grid` 위에서 ` seoul032823` 의 데이터를 내삽함 * 결과는 `seoul032823_OK`에 저장됨
head(seoul032823_OK$krige_output@coords, 2)
## LON LAT ## 1 126.50 37 ## 2 126.51 37
head(seoul032823_OK$krige_output@data$var1.pred,2)
## [1] 169.6544 169.6502
* 내삽 정보는위의 객체에 저장
myPoints <- data.frame(seoul032823)
myKorea <- data.frame(skorea)
myKrige <- data.frame(seoul032823_OK$krige_output@coords,
pred = seoul032823_OK$krige_output@data$var1.pred)
* ggplot을 사용하기 위해 데이터 형식을 data.frame으로 변형
if(!require(viridis)){install.packages("viridis"); library(viridis)}
ggplot()+ theme_minimal() +
geom_tile(data = myKrige, aes(x= LON, y= LAT, fill = pred)) +
geom_map(data= myKorea, map= myKorea, aes(map_id=id,group=group),
fill=NA, colour="black") +
coord_cartesian(xlim= LON.range ,ylim= LAT.range) +
labs(title= "PM10 Concentration in Seoul Area at South Korea",
x="Longitude", y= "Latitude")+
theme(title= element_text(hjust = 0.5,vjust = 1,face= c("bold")))+
scale_fill_viridis(option="magma")